The program is written in R.
View Code
two_arm_survival_interim_stopping = function(input_type, nsim, n, lambda1,
hratio, hratio_alt, accrual,
follow, event_frac, int_frac,
int_follow, alpha, alpha_interim,
frac_control) {
eventv1 = eventv2 = int.sev = int.eventv1 = int.eventv2 = int.betav = betav = testv = rep(NA, nsim) # storage
lambda2 = lambda1 / hratio # hazard rate for the experimental
int_accrual = min(int_frac, 1) * accrual # int_frac must be smaller or equal to 1
n1a = round(frac_control * n)
n2a = n - n1a
# If provided, convert event fraction to accrual and follow-up time
if (input_type == "Total") {
follow.v = seq(from = 0, to = follow, by = follow / 100)
accrual.v = seq(from = 0, to = accrual, by = accrual / 100)
prob1.arm1 = (accrual.v - 1 / lambda1 + (1 / lambda1) * exp(-lambda1 * accrual.v)) / accrual #arm 1 prob event
prob2.arm1 = 1 - exp(-lambda1 * follow.v) + (exp(-lambda1 * (accrual + follow.v)) / lambda1 +
(accrual - 1 / lambda1) * exp(-lambda1 * follow.v)) / accrual
prob1.arm2 = (accrual.v - 1 / lambda2 + (1 / lambda2) * exp(-lambda2 * accrual.v)) / accrual #arm 2 prob event
prob2.arm2 = 1 - exp(-lambda2 * follow.v) + (exp(-lambda2 * (accrual + follow.v)) / lambda2 +
(accrual - 1 / lambda2) * exp(-lambda2 * follow.v)) / accrual
probv.arm1 = c(prob1.arm1, prob2.arm1)
probv.arm2 = c(prob1.arm2, prob2.arm2)
probv = frac_control * probv.arm1 + (1 - frac_control) * probv.arm2 # sum events accross arms weighted by arm sample size
times = c(accrual.v, follow.v)
accstat = c(rep(1, 101), rep(2, 101))
ind1 = probv >= event_frac
indpt = min((1:202)[ind1])
if (accstat[indpt] == 1) {
int_accrual = times[indpt]
int_follow = 0
} else if (accstat[indpt] == 2) {
int_accrual = accrual
int_follow = times[indpt]
}
int_frac = min(int_accrual / accrual, 1)
}
# This is the simple case
for (jj in 1:nsim) {
# Generate full analysis data unobservables
tt1 = -log(runif(n1a)) / lambda1
tt2 = -log(runif(n2a)) / lambda2
cc1 = follow + runif(n1a) * accrual
cc2 = follow + runif(n2a) * accrual
tt = c(tt1, tt2)
cc = c(cc1, cc2)
# Generate times corresponding to the first interim
int.n1 = round(int_frac * n1a) #no pts in arm 1 at interim
int.n2 = round(int_frac * n2a)
oo1 = rev(order(cc1)) #ordering for accrual time
int.patients1 = oo1[1:int.n1]
oo2 = rev(order(cc2))
int.patients2 = oo2[1:int.n2]
int.cc = c(cc1[int.patients1], cc2[int.patients2]) - (follow - int_follow) #create true censoring times
int.tt = c(tt1[int.patients1], tt2[int.patients2]) #create true survival times
int.cc = int.cc - min(int.cc) + .01 #??? Why wouldn't I do the same for true survival times
int.status = 1 * (int.tt < int.cc)
int.group = c(rep(0, int.n1), rep(1, int.n2))
int.times = int.status * int.tt + (1 - int.status) * int.cc
int.bb = survival::coxph(survival::Surv(int.times, int.status) ~ int.group) # So we can store all the parameters we need from here
int.betav[jj] = int.bb$coef
int.sev[jj] = summary(int.bb)$coef[3]
int.eventv1[jj] = sum(int.status[int.group == 0])
int.eventv2[jj] = sum(int.status[int.group == 1])
# final analysis
status = 1 * (tt < cc)
group = c(rep(0, n1a), rep(1, n2a))
times = status * tt + (1 - status) * cc
bb = survival::coxph(survival::Surv(times, status) ~ group)
testv[jj] = bb$score
betav[jj] = bb$coef
eventv1[jj] = sum(status[group == 0])
eventv2[jj] = sum(status[group == 1])
}
# mean number of events in control/experimental arm at time of each analysis
events.int = c(mean(int.eventv1), mean(int.eventv2))
events.final = c(mean(eventv1), mean(eventv2))
# alpha = .2 # Testing Null
# alpha_interim = .05 # Testing Alternative
# Calulate Futility Stopping Probablities
# Using HR = 1 Rule
stop.betav =! (int.betav < 0)
stop.HR1 = mean(stop.betav)
# Using Test with alpha_interim
tt2 = (-((-log(hratio_alt) - int.betav) / int.sev)) < qnorm(1 - alpha_interim)
stop.testv =! tt2
stop.test = mean(stop.testv)
# Thresholds Using Test
stop.test.thresh = -(log(hratio_alt) - qnorm(1 - alpha_interim) * int.sev)
# Power with No Interim Analysis
rejectv = (testv * (betav < 0)) > qchisq(1 - alpha, 1)
power.no = mean(rejectv)
# Power with Interim Analysis - using HR1 rule
power.HR1 = mean(rejectv * (!stop.betav))
# Power with Interim Analysis - using Interim Test
power.test = mean(rejectv * (!stop.testv))
# OUTPUT
output = list(
int_frac = round(int_frac, 4),
int_accrual = round(int_accrual, 4),
int_follow = round(int_follow, 4),
n_interim_analysis = int.n1,
p_stop_int_hr = round(stop.HR1, 4),
p_stop_int_test = round(stop.test, 4),
p_stop_int_test_minus_beta = round(mean(stop.test.thresh), 4),
power_no_minus_interim = round(power.no, 4),
power_hr1_minus_interim = round(power.HR1, 4),
power_test_minus_interim = round(power.test, 4),
events_int_control = round(events.int[1], 4),
events_int_expt = round(events.int[2], 4),
events_int_total = round(sum(events.int), 4),
events_final_control = round(events.final[1], 4),
events_final_expt = round(events.final[2], 4),
events_final_total = round(sum(events.final), 4)
)
return(output)
}